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The  two-parameter  momentum-integral  method  of 
M.R.  Head  for  the  calculation  of  two-dimensional 
or  axisymmetric ,  incompressible,  laminar  boundary 
layers  is  reworked  into  an  analytical  form,  thereby 
eliminating  the  graphical  data  of  Head's  original 
method.  The  revised  method  is  coded  for  digital 
computer  solution  and  a  number  of  test  cases  are 
computed.  The  numerical  results  for  these  test 
cases  are  compared  with  very  accurate  boundary  lay¬ 
er  solutions  obtained  by  a  finite-difference  method 
and  also  results  obtained  using  the  Pohlhausen 
momentum-integral  method.  it  was  found  that  the 
Head  method  is  vastly  superior  to  the  Pohlhausen 
method  yet  requires  only  about  the  same  computation¬ 
al  expense.  Head's  method  yields  very  accurate  re¬ 
sults  even  for  highly  nonsimilar  boundary  layers. 


ADMIN  1 STRAT I VE  I NFORMAT ION 

The  work  described  In  this  report  was  carried  out  by  the  David  W. 
Taylor  Naval  Ship  Research  and  Development  Center  (DTNSRDC) ,  Ship  Perform¬ 
ance  Department,  in  support  of  the  Ship  Acoustic  Department's  Laminar  Flow 
Program.  This  program  is  sponsored  by  the  Naval  Sea  Systems  Command  (SF.A- 
037)  under  Program  Element  25634N,  Project  S0218011,  Task  Area  20052,  and 
was  performed  under  Work  Unit  1942-087. 


INTRODUCTION 

The  accurate  numerical  solution  of  the  two-dimensional  and  axisymmetr i 
incompressible  laminar  boundary-layer  equations  is  now  a  routine  matter. 
Finite  difference  methods  are  normally  used  for  such  calculations.  These 
methods  require  fairly  substantial  computing  equipment  for  efficient  imple¬ 
mentation.  Representative  examples  of  finite  difference  laminar  boundarv- 

l*  2 

layer  calculation  methods  are  given  by  Blottner  and  Cobeci  and  Bradshaw. 

Simpler  laminar  boundary  laver  calculation  methods  that  are  of  com¬ 
parable  accuracy  to  the  finite  difference  methods  are  desirable  for  several 
reasons.  For  instance,  it  would  be  useful  to  have  an  accurate  laminar 
boundary  layer  calculation  method  that  can  be  Implemented  on  a  programmable 


*A  complete  listing  of  references  is  given  on  page  45. 


desk  or  hand-held  calculator.  When  used  on  a  large  computer,  such  a  method 
would  be  considerably  less  expensive  to  execute  than  the  finite  difference 
methods . 

The  simplest  methods  of  calculating  laminar  boundary  layers  are  the 
one-parameter  momentum  integral  methods;  typically  Pohlhausen's  method  as 
described  by  Rosenhead.^  These  momentum  integral  methods  yield  acceptable 
accuracy  only  over  a  limited  range  of  variations  of  the  external  pressure 
distribution  and  are  not  always  capable  of  accurately  predicting  the  lami¬ 
nar  boundary  layer  separation  point.  Furthermore,  although  one-parameter 
momentum  integral  methods  can  be  extended  easily  to  include  extraneous  ef¬ 
fects  such  as  surface  suction  or  free-stream  unsteadiness,  the  results  com¬ 
puted  by  these  extensions  are  not  always  very  accurate. 

Two-parameter  momentum  integral  methods  for  calculating  laminar 

boundary  layers,  although  somewhat  more  complicated  than  the  one-parameter 

methods,  are,  however,  considerably  simpler  than  the  finite  difference 

methods.  Many  two-parameter  methods  have  been  considered  in  the  past  (see 
3  4 

Rosenhead  and  Waltz  for  extensive  discussions  of  momentum  integral  meth¬ 
ods).  However,  only  the  method  of  Head'  seems  to  have  fulfilled  the  expec¬ 
tation  of  considerably  higher  accuracy  in  predicting  the  boundary  layer 
properties  than  the  simpler  one-parameter  momentum  integral  methods.  Head's 
two-parameter  momentum  integral  method  (henceforth  referred  to  as  simply  the 
Head  method)  seems  to  have  received  little  attention  because  it  was  intended 
for  manual  calculation  and  thus  it  required  the  use  of  much  graphical  data. 
The  graphical  data  of  the  Head  method  is  based  on  a  seemingly  complicated 

set  of  graphical  velocity  profile  shapes  that  discouraged  analytical  de- 

3  4 

scription.  Gadd,  Jones  and  Watson  (see  Rosenhead  )  and  Waltz  give  only 
superficial  coverage  to  the  Head  method  in  their  notable  reviews  of  momen¬ 
tum  integral  methods  even  though  Head  had  shown  by  some  sample  calculations 
that  his  method  was  very  accurate. 

A  reexamination  of  the  Head  two-parameter  momentum  integral  calculation 
method  for  two-dimensional  and  axisymmetric  laminar  boundary  layers  shows 
that  his  graphical  construction  of  the  velocity  profiles  is  fairly  easily 


2 


transformed  into  an  analytical  form.  From  this  analytical  form  of  the 
velocity  profiles,  all  the  necessary  relationships  between  the  parameters 
of  the  Head  method  were  easily  computed,  thereby  replacing  ail  the  graphi¬ 
cal  data  by  polynomial  expressions.  The  method  can  then  be  easily  program¬ 
med  for  machine  calculations.  This  report  briefly  outlines  and  recasts  the 
Head  two-parameter  method  into  an  analytical  form.  The  report  contains  a 
comparison  of  some  sample  numerical  results  which  were  obtained  using  the 
Head  method,  the  Pohlhausen  method,  and  a  finite  difference  method  de¬ 
scribed  by  Blottner.^  These  results  show  the  remarkable  accuracy  of  the 
Head  method.  Although  the  modification  of  the  Head  two-parameter  method 
that  is  given  in  this  report  has  not  been  reduced  to  a  compact  form  so  that 
it  can  be  programmed  for  a  hand-held  calculator,  it  seems  that  this  could 
be  done.  Furthermore,  the  accuracy  of  the  Head  method  encourages  a  belief 
that  it  can  be  extended  for  the  accurate  calculation  of  certain  unsteady 
laminar  boundary  layers.  The  calculation  of  unsteady  boundary  layers  by 
the  two-parameter  momentum  integral  method  would  be  considerably  less  ex¬ 
pensive  than  the  calculation  of  such  boundary  layers  using  finite  difference 
methods.  A  brief  discussion  of  this  point  is  given  in  the  last  section  of 
this  report.  This  report  concludes  with  an  Appendix  in  which  the  computer 
program  for  the  Head  two-parameter  method  is  given. 

DESCRIPTION  OF  HEAD'S  METHOD 

In  this  section  a  brief  outline  of  the  Head  method  is  given.  For  the 

complete  details  of  the  method,  the  reader  is  referred  to  Head's  original 

5 

report . 

Consider  the  two-dimensional  or  axisymmetric  laminar  boundary  layer  on 

the  surface  of  an  airfoil  or  body  of  revolution  in  axial  flow.  The  problem 

is  scaled  as  follows:  Distances  are  referred  to  the  overall  length  scale  L 

of  the  body.  Velocities  are  referred  to  the  uniform  flow  velocity  at 

•y 

infinity  and  stresses  are  referred  to  the  total  head  l/2pU^  where  p  is  the 
constant  density  of  the  fluid.  Suppose  the  configuration  in  Figure  1  de¬ 
picts  a  meridian  cut  through  a  body  of  revolution  where  x  is  the  dimension¬ 
less  axial  coordinate,  s  is  the  dimensionless  surface  arc-length  starting 


Figure  1  -  Coordinates  and  Geometry  of 
a  Body  of  Revolution 


at  the  front  stagnation  point,  z  is  the  dimensionless  coordinate  normal  to 
the  body's  surface,  and  r^Cs)  describes  the  body's  profile  as  a  function  of 
the  surface  arc-length  s.  If  an  airfoil  is  being  considered  instead  of  a 
body  of  revolution,  s  and  z  can  be  interpreted  in  a  manner  similar  to  that 
shown  in  Figure  1.  The  value  of  s  is  zero  at  the  front  stagnation  point 
and  increases  downstream.  The  thickness  distribution  r^fs)  of  an  airfoil 
measures  the  distance  between  the  airfoil  offsets  and  the  camber  line. 

The  appropriate  laminar  boundary  layer  equations  can  be  found  in  any 

3 

one  of  many  standard  works  (Rosenhead,  for  example).  The  dimensionless 
laminar  boundary  layer  equations  for  planar  and  axisymmetric  external  flows 
are 
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where  U  is  the  dimensionless  Invise  Ul  slip  velocity  at  the  body's  surface 
(usually  the  potential-flow  velocity)  and  |  Is  either  0  or  l  tor  two- 
dimensional  or  axisymmettlo  t  low,  respectively.  Standard  boundary  layer 
notation  is  used  in  which  (u,  w)  denotes  the  dimensionless  tangent  tat  and 
normal  components  ot  boundary  lavei  velocities,  respectively,  p  denotes  the 
d  linens  ion  less  t  laid  pressure,  V  is  kinematic  viscosity  and  K  -  l'  l,  v  is  the 

tV 

Reynolds  number.  The  boundary  condlt ions  are 
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u  '  V  at  ~ 
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where  3  Is  the  dimensionless  boundary  layer  thickness. 

By  Integrating  Valuation  (la)  with  respect  to  /.  across  the  boundary 
laver  from  the  surface  /  0  ot  the  body  to  the  edge  ot  the  boundary  la\ei 


(s),  and  similarly  integrating  Valuation  (la)  after  first  multiplying  it  bv 


u,  one  obtains,  alter  some  algebraic  manipulations  (see  Head  or  Rosenhead  ) 
the  momentum  and  energy  Integral  equations,  respectively; 
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and  where  n  ■  z/6  and  f(n)  =  u/U.  In  deriving  Equations  (A)  it  is  assumed 
tli.it  the  velocity  profile  u/ll  is  a  function  f(n)  of  the  similarity  variable 
0  -  z/*S(s)  only,  so  that  the  partial  differential  Equations  (1)  reduce  to 
the  pair  ot  ordinary  differential  Equations  (A).  Furthermore,  the  compati¬ 
bility  cond it  ion 


t* 


di: 

ds 


m 


(b) 


is  a  consequence  of  the  requirement  that  Equation  (la)  be  satisfied  on  the 
boundary  z  0.  Thus,  the  three  Equations  (Aa) ,  (Ab) ,  and  (b)  are  used  to 


obtain  tlio  values  ol  the  two  parameters  that  determine  the  velocity  prot  lie 
t  t T> >  and  the  boundary  \aver  scale  6(s). 

The  form  of  the  velocity  profile  f(n),  that  was  given  by  Head,  is 


"  -  f(n)  -  tj(n)  +  af,(n;a)  +  bf3(n;a) 


(7) 


where  a  and  b  are  parameters,  f.(n)  1  f  g^Ol)  is  the  Hlaslus  I 'at  plate 

3 

veloeitv  profile  (Hosenhead  ) 
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Hu*  shape  functions  g.  through  g^  were  given  in  graphical  torm  bv  Head. 
These  graphs  are  not  reproduced  here.  Instead,  only  the  set  o(  analvtical 
interpolation  functions  of  these  graphs  is  given  here.  it  was  found  that 
the  functions 
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stl 


(0> 


wlu'te  p  (n) .  i  ~  1 .  r>  are  polynomials  in  n,  well  represent  t  lie  I  unc¬ 

tions  g(tn)  and  their  first  two  derivatives.  Table  1  lists  the  coetiicient 
( k) 
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hat  are  used  In  Equation  tl,l 
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TABLE  1  -  COEFFICIENTS  OF  THE  VELOCITY  PROFILE  SHAPE  FUNCTIONS 
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to. 4820091 
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561.74806 
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-305. 15578 

-434.862 
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0.0 

0.0 

-73.342991 

97. 331351 

131.01295 

With  the  aid  of  the  profile  family.  Equation  (7),  and  the  functions  of 
Equations  (8),  (9),  and  (10),  all  the  parameters  of  Equations  (5a)  through 
(5i)  can  be  evaluated  in  terms  of  the  parameters  a  and  b  as  follows.  Let 
G,  J,  E,  and  D  be  defined  by 
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where  the  upper  limit  6^/0  in  Equation  (5i)  can  be  taken  as  1  in  the 
definition  (lid) 


because  for  0^/0  >  l,  (9f/3rt)"“.  0.  Then,  by  substituting  Equations  (11a) 
through  (lid)  in  Equations  (5a)  through  (5i),  the  following  formulas  are 
obt  a  tned : 
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Kqu.it  ion  tbl  vie  Ids  the  following  formula  tor  the  variable  t*; 
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All  ot  t  hi  quantities  m  to  D*  given  hv  Equations  (12a)  through  (12e)  are 

functions  (either  ordinary  polynomials  or  rational  polynomials)  onlv  of  a 

and  b,  and  t*  is  a  function  ot  a  and  b  and  the  known  function  (dl'/ds). 

Kquat  ions  (12)  are  easily  evaluated  by  using  the  polynomials  G,  d,  E,  0. 
o  'y 

(U  e,-)^  and  (l*d/l’‘l()  which  are  listed  in  Table  2.  Thus,  Equations  (4) 
and  (hi  can  be  reduced  to  equations  determining  the  parameters  a  and  b  as  a 
function  ot  s  t ot  given  distributions  ot  the  inviscid  surface  velocity  l'(s) 
and  bodv  radius  (or  thickness)  distribution  r(1)(s). 

\t  tlu'  t ront  stagnation  point  of  blunt  bodies,  a  and  b  have  the  values 


a  ”  0.445, 


b  “  0.240  at  s  “  0 


(13a) 


anil 


a  -  0.3225,  h  -  0.1605  at  s  «  0  (13b) 
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1 


iftst 


tor  two-dimensional  and  axisymmetric  flows,  respectively.  These  values  for 
a  and  b  are  obtained  by  finding  the  roots  of  the  algebraic  equations  re¬ 
sulting  from  Equations  (4)  upon  assuming  that  both  dt*/ds  and  dh/ds  are 
bounded  in  the  limit  s  ►  0. 

2 

TABLE  2  -  FUNCTIONS  J,  G,  I),  E,  ( y- )  ,  AND  ( IN  TERMS 

*n  o  \ an L 

OF  THE  PARAMETERS  a  AND  b 


.1  -  0.32832  -  0.  31897a  +  0.14931a  -  0.1876b  -  0.14947ab 

G  -  0.12733  -  0.03003a  -  0.10089a2  +  0.13915a3  -  0.035073a'4 

+  (0.049054  +  0. IJ403a  -  0.23383a3  +  0.072643a3)b 
+  (-0.039455  +  0.093916a  -  0. 037702a2)b" 

2  3  4 

D  -  1.3713  +  1.2124a  +  l. 8754a  -  1.8268a  +  0.60524a 

+  (-0.83784  -  1.9728a  +  3. 3493a"  -  l.2919a3)b 
+  (0.94797  -  l. 5608a  +  0.6957a")b" 

K  -  0.20069  -  0.054971a  -  0. 18242a"  +  0.1758a3 
+  0.043227a*4  -  0.055451a'  +  0.0095503a6 
+  (0.060072  +  0.26478a  -  0.29465a" 

-  0.11611a3  +  0.  152a31  -  0.029892a')  b 

2 

+  (-0.  11365  +  0. 089165a  +  0. 12215a 

-  0.13673a3  +  0. 031219n4)b2  +  (0.0206b9 

-  0.049679a  +  0.040141a2  -  0. 010879a3) b3 


NUMERICAL  METHOD  FOR  SOLVING  THE  MOMENTUM  AND 
ENERGY  INTEGRAL  EQUATIONS 

In  this  section,  the  numerical  method  for  computing  the  boundary  layer 
characteristics  for  given  distributions  of  U  and  r^  is  described.  This 
numerical  method  is  implemented  in  the  computer  program  HEAD  given  in  t  he 
Appendix  and  is  Intended  for  use  on  a  digital  computer. 

The  numerical  solution  of  Equations  (4)  can  be  approached  in  two  wavs. 
The  first  way  is  by  the  direct  calculation  of  t*  and  h  using  Equations  (4). 
This  would  require  the  solution  of  a  system  of  algebraic  equations  for  the 
parameters  a  and  b  in  terms  of  the  known  values  of  t*  and  h  at  each  step 
of  the  calculations.  The  second  way  to  numerically  solve  Equations  (4)  is 
to  rewrite  them  in  terms  of  the  parameters  a  and  b.  This  eliminates  the 
necessity  of  solving  algebraic  equations  at  each  step  and  is  the  basis  of 
the  numerical  procedure  described  below. 

Equations  (4)  can  be  written  in  the  form 


da 

ds 


F(.i.b) 


db 

ds 


0( a, bl 


where 


f  (a ,  b) 


G(a , bl 


A 


mU”  2 
U'  U' 


mU(r  V 

l'1  + 

r  ^ 

0 


-U*  { 2D*-h [  t-Hn (H+ 1 ) ) } 
mU 


V 


2m  3h  _  2m  2h 
'  3a  3b  3b  3a 


(14a) 


(lsb) 


i  1  Sal 


(15b) 


( l  Sc  2 


( 1  Sd  2 


1 1  Sol 


1 1 


w 


and  where 


0*  f 

ds 


and , 


<r0J) 


dr 


J 


0 

ds 


(15f) 


Equations  (14)  can  be  integrated  numerically  using  the  predictor- 
corrector  method  that  is  based  on  the  generalized  midpoint  and  trapezoidal 
quadrature  formulas.  The  values  of  a(s)  and  b(s)  at  the  points  s^  of  a 


general  grid  {s.:  I  »  0,  1,  2,  ...  are  denoted  by  a^  and  b^,  respectively. 
Then,  given  the  values  a^  and  b^  of  a  and  b,  respectively,  at  s^,  the 
values  a.^  and  b,+j  at  sj+j  are  obtained  using  the  formulas 


i  +  1 


l+l 


ai+l  '  ai 


bi+l  *  bi 


i-1 

*  i—  1 

As 

+  — 

As 


‘i  “”t' 


/  ,  '  .  in  ,  ' 

i  t 


i  +  1 


1+1 


i> 

(16a) 

i> 

(16b) 

)) 

(17a) 

)) 

(17b) 

where  As^  +  j  •  *>j+|  -  s^  and  c  *  As.^/As.  for  i  m  2,  3,  ...  For  the  first 
station  s^  past  the  stagnation  point  *  0,  a  special  procedure  is  used. 


In  order  to  avoid  the  front  stagnation  point  where  r^  = 


0  and  da/ds 

and  db/ds  are  indeterminate,  the  numerical  integration  is  started  by 
Iterating  the  backward  Euler  formulas  as  follows: 


(k)  x  a  «r/0(k-l)  K(k_1>  o  N 

a^  -  aQ  +  As^  r (a j  ,  bj  ,  s^ 


(18a) 


* 


(k)  ,  i  a  n,  (k~D  u(k-D  \ 

>j  *  by  +  Asj  G(a|  ,  bj  ,  s^) 


(18b) 


where  k  -  1,  2,  3,  ....  aj°^  -  aQ,  b|0)  -  b0  and  the  third  or  fourth 
(k  -  3  or  4)  iterates  usually  are  sufficiently  accurate  approximations  of 
a^  and  b^. 


i. 
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The  numerical  integration  method.  Equations  (16)  and  (17)  changes 
the  step  size  As^+^  as  the  integration  proceeds  downstream  from  the  point 
s  =  Sp  thus  automatically  generating  the  grid  {s^}  in  such  a  way  as  to 
keep  a  certain  measure  e  of  the  numerical  Integration  error  below  a  preset 
value  E.  Each  new  step  s^+^  =  s^  +  As^+^  of  the  numerical  integration  pro¬ 
cedure  is  obtained  from  the  previous  step  =  s^  ^  +  As^  by  setting  As^+.j  = 
AAs^.  In  the  first  attempt  at  moving  from  step  s^  to  s^+^,  A  is  either  set 
equal  to  1  or  is  the  value  given  by  the  procedure,  described  below,  for 
increasing  the  step  size  if  higher  accuracy  than  desired  is  being  obtained. 

Suppose  the  values  of  a.,,,  b.,,,  a.,.,  b,  ,  have  been  computed  at 
‘  l+l  l+l  l+l  i+1 

s.,,  =  s.  +  AAs . .  Let  e  be  the  measure  of  the  absolute  error,  defined  by 
l+l  i  i 


la 


e  =  max 


i+1  l+l 


1  bi+l~bi+l I  1 


(19) 


where  d  =  3  +  2/c^.  If  e  >  E  the  current  values  of  s..,,  a..,,  a..,, 

l+l  l+l  l+l 

and  b^+^  are  discarded  and  a  new  value  of  A  is  estimated  by  using  the 
formula 


1/3 

) 


b. 


i+1 


(20) 


and  new  values  of  s .  ,  ,  =  s.  +AAs.,  a.,.,  a.,,,  b.,,  and  b  ,  are  computed . 
i+l  i  l  l+l  i+1  l+l  l+l 

The  procedure  is  repeated  until  the  inequality  e  <  E  is  satisfied.  In  all 
the  numerical  experiments  using  this  procedure  to  satisfy  the  accuracy  re¬ 
quirement  that  e  ^  E  it  has  never  been  necessary  to  repeat  the  procedure 
more  than  once  at  each  step. 

If  too  much  accuracy,  say,  for  example,  e  v  E/10,  is  being  obtained 

at  step  s,,,  then  the  current  value  of  As  is  accepted  but  the  next 
i+1  i+l 

location  s.+  ,  is  obtained  by  using  the  value  As^+,  =  AAs^+j  where 


(21) 
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This  method  of  changing  the  integration  step  size  has  proven  to  be  very 
effective  in  many  numerical  experiments.  The  value  of  L  that  has  been  used 
in  all  the  computations  presented  in  this  report  is  E  »  10  J. 

Special  care  must  be  exercised  when  U'  vanishes.  Founula  (12f)  shows 
that  when  U'  “  0,  m  must  also  be  zero  in  order  for  t*  to  be  finite.  When¬ 
ever  U '  is  close  to  zero  (say  %  10  )  at  location  s^+^,  a  special  procedure 
is  used  to  evaluate  t*  =  -tn/U'  in  order  to  maintain  accuracy.  Instead  of 
simply  evaluating  the  ratio  m/U'  in  order  to  obtain  t*  at  the  point  s  , 
t*i+1  is  obtained  by  numerically  integrating  Equation  (4a)  vising  a  special 
explicit  procedure  based  on  the  known  values  t*._j,  t*^  (and  the  rest  of 
the  data  needed  in  Equation  (20))  at  the  points  s.  ^  and  s^,  respectively. 
The  value  t* of  dt*/ds  at  the  point  s.  can  be  computed  using  Equation  (4a) 
Then,  the  value  of  t*^+^  at  the  point  s can  be  computed  using  the  two- 
point  explicit  Hermlte  extrapolation  formula 


i+1 


‘Vl  + 


As. 


-  ‘*1 


(ASi+As1+1) 


As. 


(r  *  _t*  ) 

1  i-i  i-r 

o 

As  ~ 


(Asi+As1+i)- 


(22) 


If  tlie  values  of  t*.  t*j  and  the  values  of  the  rest  of  the  data  needed 

in  Equation  (4a)  are  known  to  second-order  accuracy  in  the  step  sizes  As^, 

then  the  value  of  t*...  computed  using  Equation  (22)  is  of  second-order 

l**’  l 

accuracy  in  the  step  sizes  As,  and  As .  ,  ,  .  Thus,  even  when  1!'  and  m  become 

i  i+1 

very  small,  their  ratio  t*  can  still  be  calculated  accurately  using  Equation 

(22). 


The  characteristics  of  the  boundary  layer  are  determined  once  the  dis¬ 
tributions  along  the  body's  profile  of  the  values  of  a  and  b  are  known. 

The  boundary  layer  parameters  that  are  usually  of  main  interest  are  the  dis¬ 
placement  thickness  =  6J,  the  momentum  thickness  8  =  AG,  the  shape  fac¬ 
tor  H  »  6./0,  and  the  skin  friction  coefficient  cf  =  T  /  [4(pU„  )1  where  r 

1  I  w 


is  the  dimensional  wall  shear  stress  pv(9u/8z)q.  These  parameters,  scaled 
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m  ■■ 


only  by  the  overall  body  length  L,  overall  velocity  Uq  and  the  fluid  density 
P  are  given  in  terms  of  the  calculated  values  of  t*  and  H  by 


(23a) 


6  j  *  0H 


(23b) 


2U£ 


(23c) 


The  intrinsic  boundary  layer  scale  5(s)  can  be  computed  using  the  relation 


once  the  distributions  of  a  and  b  are  known.  Thus,  the  boundary  layer 
velocity  profiles  f (n)  “  u/U  can  be  obtained  using  Equations  (8)  through 
(10).  The  velocity  profiles  f(n)  can  be  plotted  in  terms  of  any  other  body 
normal  coordinate  n.  that  is  scaled  by  some  other  length  scale  L,  by  simply 
multiplying  n  by  Li5/L  to  obtain  n. 


NUMERICAL  RESULTS  AND  DISCUSSION 

The  laminar  boundary  layer  characteristics  were  calculated  for  three 
sample  bodies  of  revolution  by  three  methods;  the  Pohlhausen  and  Head  mo¬ 
mentum  integral  methods  and  the  finite-difference  method  designated  DB, 
that  is  described  and  referred  to  as  the  Davis-B  scheme  by  Blottnor.  The 
DB  finite-difference  method  provides  second-order  accuracy  in  the  step 
sizes  As  and  An  along  and  normal  to  the  body,  respectively.  The  values  of 
the  step  sizes  As  and  An  that  were  used  in  the  DB  method  were  small  enough 
for  each  body  considered  in  this  report  that  further  reductions  in  the  step 
sizes  made  less  than  a  one-percent  difference  in  any  of  the  computed  values 
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of  the  boundary  layer  characteristics.  Accordingly,  the  results  obtained 
using  the  OB  method  are  presumed  to  be  exact  for  comparison  purposes.  As 
a  further  check  of  these  computational  results,  the  boundary  layer  char¬ 
acteristics  were  also  calculated  for  several  bodies  using  a  computer  pro- 
gram  based  on  the  Kel ler-Cebec i  box  method."  The  numerical  results  obtained 
with  the  box  method  in  all  cases  agreed  to  within  two  percent  of  the  results 
obtained  with  the  OB  method. 

The  boundary  layer  characteristics  were  computed  by  all  three  methods 
tor  a  number  of  bodies  (both  cylindrical  and  axisymmetric  bodies).  It  is 
sufficient  tor  comparison  of  the  accuracy  and  computational  cost  of  the 
Pohlhausen,  Head,  and  DB  methods  to  illustrate  the  numerical  results  on 
only  three  sample  bodies  since,  in  all  the  cases  considered.  Head's  method 
gave  equallv  accurate  results.  The  three  illustrative  bodies  that  were 
chosen  are  fairly  extreme  cases  for  which  the  prediction  ol  the  laminar 
boundary  layer  can  be  classified  as  either  fairly  easy  or  hard.  The  easiest 
possible  case  is  probably  the  flat  plate  boundary  layer  but  this  was  not 
deemed  a  fair  test  because  Head's  method  has  an  accurate  interpolation  of 
the  flat  plate  laminar  boundary  layer  velocity  profile  built  into  it. 

The  two  easy  cases  considered  for  this  report  are  the  flow  past  a 
sphere  and  the  axial  t low  past  a  spheroid  with  a  length-to-d iameter  ratio 
of  4  to  L.  These  are  common  test  cases  for  which  many  approximate  boundary 
layer  solutions  have  been  published.  Any  viable  calculation  method  must  be 
able  to  predict  the  laminar  boundary  layer  on  these  bodies  with  reasonable 
accuracy. 

The  hard  test  case  considered  for  this  report  is  a  rather  contrived 
mathematically  defined  body  of  revolution  with  a  flat  nose.  Such  extremely 
blunt  bodies  have  very  rapid  and  large  variations  of  the  surface  pressure 
distribution.  These  variations  can  cause  approximate  laminar  boundary  layer 
calculation  methods  to  predict  catastrophically  inaccurate  boundary  layer 
characteristics  such  as  separation  where  it  does  not  occur  or  no  separation 
where  it  does  occur.  The  hard  test  case  chosen  for  illustration  is  desig¬ 
nated  body  B.  On  body  B,  the  boundary  layer  nearly  separates  at  a  point 


near  the  nose  of  the  body  but  recovers  and  eventually  culminates  in  actual 
boundary  layer  separation  near  the  tail  of  the  body.  It  seems  that  body  B 
is  a  sufficiently  extreme  case  for  laminar  boundary  layer  calculation  meth¬ 
ods  so  that  if  a  method  is  successful  for  body  B,  then  one  can  expect  it  to 
be  successful  on  most  other  bodies  (two-dimensional  or  axisymmetric  ones). 

The  first  test  case  to  be  considered  is  the  boundary  layer  flow  on  a 
sphere.  It  is  sufficient  to  compare  only  the  surface  distribution  of  the 
values  of  the  shape  factor  H  and  skin  friction  coefficient  C^Vr^.  Figures 

2  and  3  show  the  distributions  of  the  surface  values  of  and  of  H, 

respectively,  plotted  as  a  function  of  the  axial  coordinate  x.  The  front 
stagnation  point  is  located  at  the  point  x  =  -1.  The  values  of  and 

H  that  are  obtained  by  the  Pohlhausen,  Head,  and  DB  methods  are  denoted  by 
the  triangular,  circular,  and  diamond  shaped  symbols,  respectively.  This 
notation  is  used  in  all  the  comparisons  shown  in  this  report.  On  the  front 
of  the  sphere,  from  the  stagnation  point  at  x  =  -1.0  to  the  location  x  = 
-0.3,  the  numerical  results  that  were  obtained  by  the  Pohlhausen,  Head, 
and  DB  methods  are  plotted  at  the  different  station  values  x  for  clarity. 
From  the  station  x  =  -0.3  to  the  point  of  separation  near  x  *  0.26  the 
numerical  results  for  and  H  obtained  by  each  of  the  three  methods  are 

plotted  at  the  same  station  values  x. 

It  can  be  seen  in  Figures  2  and  3  that  the  laminar  boundary  layer  cal¬ 
culation  method  of  Head  is  considerably  more  accurate  than  the  Pohlhausen 
method  although  the  latter  seems  to  be  adequate  for  rough  estimates.  The 
separation  point  is  predicted  by  the  DB  method  to  occur  at  the  location 
x  =  +0.255.  Head's  method  predicts  the  separation  point  at  x  =  0.265  and 
Pohlhausen's  method  predicts  it  at  x  =  0.305.  The  price  that  one  has  to 
pay  for  Head's  method,  over  the  price  of  the  Pohlhausen  method,  for  the 
more  accurate  boundary  layer  predictions  is  very  low.  The  calculations  for 
the  sphere,  given  by  the  computer  program  of  the  Head  method  in  the  Appen¬ 
dix,  required  3.5  seconds  to  execute  on  the  CDC  6700  computer.  The 
Pohlhausen  method  required  4.7  seconds  using  a  less  efficient  integration 
routine.  A  comparable  numerical  integration  method  to  the  one  used  for 
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Hoad's  mot  hod  could  probably  reduce*  by  halt  the  execution  time  of 
1’oh  1  hausen '  s  method.  Thus,  there  is  little  to  gain  on  a  computer  by  using 
Poh 1  hausen 1 s  method  in  favor  of  Head's  method.  Furthermore,  there  is  little 
ditference  in  computer  storage  requirements  between  Pohl hausen ' s  and  Head's 
method,  but  this  will  be  discussed  later.  The  OH  method  required  25  seconds 
to  execute  for  the  sphere.  (The  Cebeci-Kel ler  box  method  required  slightly 
more  execution  time  than  the  DB  method  for  all  the  bodies  cons idered . ) 

l'he  comparison  of  the  surface  distributions  of  values  of  C  » /R^  and  II 
on  the  spheroid  as  calculated  by  the  Pohlhausen,  Head,  and  1)B  methods  are 
shown  in  Figures  4  and  5.  Kssentially  the  same  comments  can  be  made  con¬ 
cerning  these  comparisons  as  the  ones  for  tin*  sphere.  The  Pohlhausen  pre¬ 
dictions  for  the  spheroid  are  better  than  for  the  sphere,  although  still 
not  as  good  as  Head's  predictions.  line  must  not  be  misled  by  the  seemingly 
better  predictions  of  the  values  of  11  near  separation  on  the  spheroid  by 
Pohlhuusen's  method  than  by  Head's  method.  It  should  be  noted  that  the 
graph  of  11  predicted  by  Pohl hausen ' s  method  crosses  the  exact  DB  graph  ol 
H  very  close  to  separation  and,  thus,  the  good  agreement  between  values  ot 
11  there  is  somewhat  tortuitous.  In  contrast,  separation  (near  point  x 
O.bbS  in  Figure  4)  is  predicted  slightly  more  accurately  by  Head's  method 
than  by  Poh 1  hausen ' s  method.  The  computer  execution  times  for  calculating 
the  boundary  layer  characteristics  on  the  spheroid  were  1,  7,  and  IS  seconds 
for  the  Head,  Pohlhausen,  and  DB  methods,  respectively. 

The  final  test  case  tor  comparing  the  Laminar  boundary  layer  calcula 
t ion  methods  is  body  B  which  is  described  by  the  following  set  ot  oqoat  Ions. 
The  flat  nose  is  given  by 

x  -  0,  for  re  [0.0,  0.025]  (25a) 


the  spheroidal  blending  region  is  given  by 


0.025  +  0.021)7 


,  1/2 


'.111  (0.II3) 


(2  5b) 


for  xtiO.O,  0.113] 
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Computed  Skin  Friction  Coefficient  on  a  ^  to  1  Spheroid 


tho  para  I  lot  middle  body  Is  given  by 


»■  “  0.04637  (2'h-) 

for  xt  [0.113,  0.592) 

and  tho  tall  piece  is  given  by 

r()  -  0.09274  (z.  +  1.137153Z"2  -  10.774885z'  +  19.78428bzA 

-16.  792534*  1  +  5.645977zh) 1  2  (25d) 


for  xt  [0.592,  1.01,  7-  ”  1.47049x  -  0.47049 

ll  should  bo  notod  that  this  prot  i  1 «'  equal  ion  lor  body  8  is  used  only 
tv>  compute  a  sot  ot  discrete  points  { x ^ ,  i^fx^)}  describing  tho  hodv.  This 
sot  ot  discrete  points  Is  used  as  input  to  calculate  tho  hodv  surface  pot  on 
tial  t  low  velocity  values  i  l  ’  ^  }  at  a  sot  ot  points  midwav  between  the  dot  til¬ 
ing  points  Ir^tXj)}.  lhe  values  ot  the  potential  t  low  velocities  {l^!  are 

(l 

computed  using  a  pane  l  method  developed  at  ONTSKDC  by  t'hang  and  I’ ten.  the 
laminar  boundary  lavei  programs  I  ot  the  I’ohlhuusen,  Head,  and  Oft  methods 
that  were  used  for  the  calculations  described  in  this  report,  used,  as  in¬ 
put,  the  art  of  values  (r^(Xj)l  and  il'^)  and  derived  cubic  spline  interpo¬ 
lation  potynominals  t  or  the  distributions  ol  r^(s)  and  U(s).  1‘hiis,  tin- 
actual  body  B  used  t or  the  csleulat  ions  is  the  cubic  spline  t  it  to  the  data 
points  (r^fx.))  that  l«  implicit  in  the  boundary  Inver  programs.  The  cubic 
spline  interpolat  ion  polynomials  I  oi  i  ^tsl  and  Ilfs)  at  <-  exactly  the  same 
I  or  all  three  boundary  layer  programs. 

I'be  distribution  ot  values  ot  the  surlaee  pressure  coot  t  (clout  t'  on 
hodv  B  is  plotted  versus  surlaee  arc  length  s  in  Figure  6.  i'be  trout  stag 


nation  point  corresponds  t  o  the  point  s  if.  0 .  Note  the  extremely  steep 

slope  ol  the  graph  ot  F  near  the  nose  amt  its  sudden  reversal  neat  the 
point  s  0.01.  This  rapidlv  changing  pressure  distribution  near  the  nose 
ot  hodv  B  can  be  expected  t o  cause  large  dovlat ions  ot  the  local  bouudarv 
1  aver  characteristics  from  the  cliarae  t  or  1st  ies  of  stmllaiitv  boundary  lavei 


flow.  Thus,  one  would  expect  most  low  order  momentum  integral  methods  to 
have  difficulties  in  accurately  predicting  the  boundary  layer  character¬ 
istics  for  body  B. 

Figures  7  and  8  show  the  remarkably  accurate  prediction  of  the  values 
of  and  H  by  Head's  method.  The  values  of  C^/Rj  and  H  near  the  nose 

iO.O  s  v  0.  1  for  and  0.0  <  s  <  0.2  for  H)  were  plotted  separately 

from  the  rest  of  the  body  in  Figures  7a  and  8a,  respectively,  because  of 
their  very  large  magnitudes  and  rapid  variations  there.  The  accuracy  of 
the  predictions  of  the  values  of  Cj.*K  and  H  by  Head's  method  are  even  more 
remarkable  in  the  light  of  the  complete  failure  on  body  B  of  Pohlhausen's 
method.  The  predicted  values  of  the  Pohlhausen  parameter  A  were  outside 
the  region  in  which  the  Pohlhausen  velocity  profiles  are  meaningful  approx¬ 
imations  of  the  exact  velocity  profiles.  Therefore,  the  Pohlhausen  method 
broke  down  and  could  not  be  continued  past  the  point  s  4  0.02f>.  The  same 

type  ot  breakdown  of  the  Pohlhausen  method  occurred  at  the  same  point, 

s  0.02t>,  using  two  different  types  of  numerical  Integration  methods  on  a 
number  of  different  fixed  and  variable  integration  nets  Is.}.  Accordingly, 

l 

it  was  concluded  that  the  Pohlhausen  method  is  incapable  ol  reasonably 
approximating  the  boundary  layer  on  the  flat  lace  of  the  body  B.  It  is 
interesting  to  note  how  accurately  Head's  method  predicts  the  laminar 
separation  point  In  Figure  7b,  especially  in  view  of  the  severe  pressure 
variations  that  the  boundary  layer  had  to  negotiate  before  arriving  at  the 
actual  laminar  separation  point.  In  particular,  it  is  to  be  noted  that  t he 
laminar  boundary  layer  nearly  separates  at  the  point  s  -  0.145  on  body  B 
but  then  rapidly  recovers  downstream  of  this  point.  Head’s  method  faith¬ 
fully  reproduces  this  behavior  of  the  boundary  layer  on  body  B.  The  rela¬ 
tive  errors  of  the  values  of  (1  /r^  obtained  by  Head's  method,  i.e., 

((V'Vhead  -  (V‘Vi)B)/(Cf^VnB  aro  f;,lrly  snu,U’  le88^_um  s  pcrcent 

except  near  some  Isolated  points  where  the  slope  ot  the  i'(tRj  versus  s 
curve  Is  nearly  vertical  at  the  point  where  the  boundary  layer  nearly  sepa¬ 
rates.  Although  the  relative  errors  at  the  latter  points  are  large,  it  is 
not  deemed  serious  because  the  absolute  value  of  ('  j- »  K  ^  is  small  there. 
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Figure  8  -  Computed  Shape  Factor  Parameter  on  Body  B 


It  is  interesting  tu  note  that  for  the  case  of  body  B,  Head's  method 
required  19  seconds  of  execution  time  and  the  DB  method  required  217  sec¬ 
onds.  By  requiring  less  accuracy  of  the  DB  method,  for  example,  accuracy 
comparable  to  the  results  of  Head's  method,  the  step  sizes  required  by  the 
DB  method  can  be  increased  to  large  values  that  require  only  about  half  the 
execution  time  of  the  accurate  step-size  DB  method.  Thus,  for  equivalent 
accuracy  levels,  Head's  method  would  still  be  about  five  times  as  fast  as 
the  DB  method.  It  is  also  interesting  to  note  that  the  automatic  step  size 
selection  mechanism  of  the  numerical  method  that  was  used  to  implement 
Head's  method  reduced  the  step  sizes  As.  to  values  on  the  order  of  10  5  in 

t 

order  to  accurately  integrate  the  boundary  layer  Equations  (14)  near  the 

nose  of  body  B.  On  the  parallel  middle  body  portion  of  body  B,  the  step 

_2 

sizes  As.  automatically  increased  to  values  on  the  order  of  10  . 

l 

There  are  several  reasons  why  it  may  also  be  necessary  to  compute  the 
boundary  layer  velocity  profile  f (q)  =  u/U  and  its  first  two  derivatives, 
e.g.,  for  boundary-layer  instability  analysis.  Thus,  a  comparison  of  the 
computed  values  of  the  velocity  profile  f  and  its  first  two  derivatives  f' 
and  f"  that  were  obtained  from  Head's  method  and  the  DB  method  is  made  in 
Figures  9  through  11.  These  velocity  profiles  are  plotted  versus  the  co¬ 
ordinate  n*  that  is  normal  to  the  body  surface  and  that  is  made  dimens ion¬ 
less  by  the  local  value  of  the  displacement  thickness  6.  (s).  Figures  9 

through  11  show  the  body  B  velocity  profiles  and  their  derivatives  df/drf* 

2  2 

and  d^f/drij.^  at  locations:  s  *  0.026,  where  C  is  a  minimum  and  is  rapidly 

*  p 

varying;  at  s  '  0.14b,  where  boundary  layer  nearly  separates;  and  at  s 
0.862,  where  the  boundary  layer  is  just  beginning  to  separate.  The  values 
of  the  velocity  profiles  f (n)  are  very  accurately  predicted  by  Head's  meth¬ 
od  at  all  of  the  above  three  points  s  on  the  body.  The  values  of  the  first 
derivatives  df/dru  are  less  accurate  but  still  fairly  good.  Substantial 
inaccuracy  is  found  only  in  the  values  of  the  second  derivative  d~f  dn*  . 
The  major  discrepancies  between  the  values  of  d  f/dr)ft~  predicted  by  Head  s 
method  and  the  DB  method  are  mainly  confined  to  a  region  very  close  to  the 
body  surface. 
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Figure  LOh  -  Comparison  of  I  1  (n*)  and  i"(n*)  1'rotiles 


Figure  10  -  Computed  Velocity  Profiles  on  Rody  R  at  s 
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It  is  really  remarkable  how  accurately  Head's  method  predicts  the 

velocity  profiles  f(q).  It  is  probably  impossible  to  obtain  much  better 

predictions  of  a  general  velocity  profile  shape  f (n)  and  its  derivatives 

bv  any  other  two-parameter  method.  This  can  be  surmized  by  noting  in 

2  2 

Figures  10b  and  11b  that  the  second  derivative  d~f/dn*  has  several  local 
minimums  and  maximum,  across  the  boundary  layer.  Such  variable  functions 
would  seem,  in  general,  to  require  more  than  two  parameters  to  accurately 
interpolate. 

CONCLUDINC.  REMARKS 

The  numerical  results  of  the  laminar  boundary  layer  calculations  using 
the  Pohl hausen.  Head,  and  DB  methods  shown  in  Figures  2  through  11  show 
that  Head's  method  is  very  accurate  and  vastly  superior  to  Pohlhausen's 
method.  Since  Head's  method  requires  no  more  computer  execution  time  and 
very  little  more  computer  core  storage  than  Pohlhausen's  method,  there 
seems  little  reason  for  using  the  latter  method  when  the  former  is  avail¬ 
able.  However,  with  the  availability  of  accurate  finite-difference  meth¬ 
ods,  such  as  the  DB  method,  what  advantages  are  there  in  using  Head's 
method?  If,  for  example,  extremely  high  accuracy,  essentiaLlv  exact,  lami¬ 
nar  boundary  layer  computations  are  required,  then  one  must  use  the  finite 
difference  methods.  However,  if  only  moderately  accurate  results,  as  ob¬ 
tained  from  Head's  method  shown  here,  are  required,  then  Head's  method  has 
several  advantages. 

The  first  advantage  is,  of  course,  that  Head's  method  requires  onl\ 
about  one-fifth  the  computer  execution  time  of  the  DB  and  similar  1  Lnite- 
difference  methods  for  the  same  level  of  accuracy.  Although  computer  exe¬ 
cution  times  are  not  large  even  for  the  finite-difference  methods,  the 
savings  gained  by  using  Head’s  method  can  be  substantial  for  design  where  a 
parametric  study  of  the  boundary  layer  behavior  on  a  large  number  of  bodies 
must  be  made. 

Secondly,  there  is  a  substantial  difference  in  computer-core  storage 
requirements  between  Head's  method  and  finite-difference  methods.  In  fact, 


11 
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the  small  computer-core  storage  requirements  of  Head's  method  raise  the 
possibility  that  with  further  refinements  the  method  can  probably  be  exe¬ 
cuted  on  a  programmable  hand  calculator.  However,  for  large  computers, 
there  is  also  an  advantage  of  the  small  computer-core  storage  requirements 
of  Head's  method  when  one  needs  complete  information  of  the  boundary  layer 
characteristics  over  the  whole  body,  including  boundary  layer  velocity  pro¬ 
files  at  each  of  many  locations  on  the  body.  For  finite-difference  methods 
such  a  requirement  would  necessitate  the  storage  of  many  vector  arrays  that 
describe  the  velocity  profiles  by  points.  But  Head's  method  requires  only 
the  storage  of  the  two  arrays  {a^}  and  {b^}  of  the  values  of  the  parameters 
a  and  b  at  the  locations  {s^}  i  =  0,  1,  2,  3,... I  on  the  body.  Furthermore, 
if  the  boundary  layer  characteristics  for  a  certain  body  must  be  stored 
permanently  for  some  future  purpose  it  is  very  easy  to  simply  punch  the 

arrays  {a.}  and  {b.}  on  a  thin  card-deck, 
l  i 

Since  it  has  been  found  that  Head's  method  is  very  accurate,  it  seems 
that  the  method  can  be  extended  to  accurately  predict  unsteady  boundary 
layers  of  low  and  moderate  values  of  frequency  (or  low  and  moderate  values 
of  inverse  time  scale) .  Such  an  investigation  is  worthwhile  and  should  be 
fairly  easy  to  carry  out.  The  unsteady  versions  of  Equations  (14)  can  be 
shown  to  have  the  form 


.  9a  _9a  _3b  ,  _  lil 

A11  3t  +  A12  3s  +  B11  3t  B12  3s 


=  C, 
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3s 


=  C„ 


(26a) 

(26b) 


where  the  coefficients  A..,  B..  and  C.  are  functions  of  a,  b,  r_(s)  and 
U(s,t).  Equations  (26)  are  in  a  natural  form  to  which  one  can  apply  the 
method  of  characteristics.  The  method  of  characteristics  is  probably  the 
most  effective  available  numerical  method  for  solving  a  pair  of  first  order 
hyperbolic  equations  in  two  independent  variables.  In  contrast  to  this,  a 
finite  difference  time  dependent  calculation  of  the  momentum  Equation  (la) 
involves  three  independent  variables  and  thus  would  require  much  larger 
storage  arrays  and  much  more  computer  time. 
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APPENDIX 

THE  COMPUTER  PROGRAM 


The  FORTRAN  IV  computer  program  that  Implements  Head's  two-parametet 

two-dimensional  and  axisymmetr ic ,  incompressible  laminar  boundary  laver 

calculation  method  is  listed  in  this  appendix.  The  program,  called  HEAD, 

requires,  as  input  for  an  axisymmetric  body,  a  discrete  set  of  offsets 

ix.,r,,  }  where  {x.}  is  the  axial  coordinate  and  {r,.  }  is  the  radius,  and 
i  0 ,  i  0 , 

i  i 

the  invisc id  slip  velocity  lU,}.  The  arc-length  {s.}  is  computed  by  the 
program.  For  a  two-dimensional  body,  the  sets  of  arc-length  values  {s.  i 
and  invisc id  slip  velocities  (U. !  are  the  required  input.  Also  for  two- 
dimensional  cases,  the  arc-length  location  s(1  of  the  front  stagnation 
point  must  be  provided. 

The  program  HEAD  is  presently  set  up  to  run  on  the  Burroughs  B7000 
computer  directly  in  the  interactive  mode  from  a  terminal.  Thus  all  one 
has  to  do  to  run  the  program  is  to  give  the  command  to  EXECUTE  HEAD.  The 
input  must  be  typed  in  at  the  terminal,  in  free  format,  in  the  order  in 
which  the  program  calls  for  it.  In  each  case,  the  program  prompts  the 
user  as  to  the  input  that  must  be  provided.  For  example:  the  first  input 
card  required  by  the  program  is  the  set  of  parameters  NT.  IPOTF.  SCAM', 
IAN.  flic  program  first  types  a  message  requesting  that  these  parameters 
be  defined.  NT  is  the  number  of  input  offsets  defining  the  body 
(1  <  NT  <  200);  IPOTF  equals  zero  (0)  for  axisymmetr ic  flow  and  one  (11 
for  two-dimensional  flow;  SCALE  is  an  arbitrary  scale  factor  to  rescale 
the  input  offsets  at  the  discretion  of  the  user  (normally  SCALE  =  1.0), 

IAN  equals  zero  (0)  if  the  body  and  potential  flow  is  provided  by  t  lie 
user,  equals  one  (1)  if  the  user  wishes  to  use  a  built-in  potential  flow 
for  the  sphere,  and  equals  two  (2)  it  the  user  wishes  to  use  tin'  built-in 
potential  flow  for  a  circular  cylinder  with  no  circulation  (the  front 
stagnation  point  is  at  0.0).  In  the  latter  two  cases  one  can  set  NT— 1 
and  input  for  body  offsets  and  potential  flow  velocities  are  not  required. 


The  Input  variable  "initial  step  size"  for  which  the  program  will  ask 
serves  a  dual  role.  The  first  role  is  precisely  that  of  setting  the 
first  integration  step  at  the  front  stagnation  point.  However,  this  is 
not  important  and  the  program  is  not  sensitive  to  the  value  of  initial 
step  size  since  it  will  automatically  be  adjusted  to  maintain  a  certain 
integration  accuracy.  The  second,  more  important,  role  of  initial  step 
size  is  that  it  determines  the  approximate  interval  along  the  body  surface 
at  which  output  is  provided.  Near  separation,  more  output  may  be  provided 
than  far  in  front  of  separation.  Thus,  the  density  of  the  output  of 
boundary  layer  characteristics  along  the  surface  of  the  body  is  controlled 
by  the  number  put  in  under  the  heading  "initial  step  size." 

The  HEAD  program  yields  the  following  output.  First,  a  summary  of 
the  geometric  (r^(s))  and  inviscid  flow  (U(s))  data  is  printed  under  the 
columns  of  R,  S,  UB,  A,  R1  which  are  the  input  values  of  body  radius  r^, 
the  axial  locations  of  the  input  values  of  body  radius,  the  input  values 
of  U,  the  arc-length  values  of  s  at  which  the  input  values  of  U  are  speci¬ 
fied,  and  the  body  radius  r ^  at  the  positions  of  the  input  values  of  U, 
respectively.  In  the  version  of  the  HEAD  program  given  in  this  report, 
the  input  values  of  U  correspond  to  locations  along  the  body  which  are 
midway  between  consecutive  locations  at  which  the  input  values  of  r ^  are 
given.  Following  the  above  output,  the  boundary  layer  characteristics  are 
printed.  The  numerical  values  of  the  quantities  arc-length  location  s . , 
step  size  As  ,  body  radius  r  ,  first  derivative  with  respect  to  s  of  U  at 

s.,  second  derivative  with  respect  to  s  of  U  at  s.,  H.,  (Cr»R, ).,  (o,vR, ),, 
j  jj  fLj  lLj 


and  (OvRj).  are  printed  under  the  column  headings  x,  HI,  ROB,  UB,  C,  D, 
DUBDS,  DUBDS2,  H,  CT,  DELST,  and  J,  respectively.  The  program  terminates 
either  at  the  end  of  the  body  or  just  in  front  of  separation  with  the 
message  that  the  step  size  of  the  numerical  integration  procedure  is  too 
small  (0(10~7)). 


In  this  program  no  provisions  for  computing  the  boundary  layer 
velocity  profiles  have  been  made.  Given  values  of  a  and  b  (t.e.,  C  and  D 
in  the  program)  one  can  easily  write  programs  using  formulas  (7)-(10)  to 
construct  velocity  profiles. 
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: 


290  C*******«PROGRAN  HEAD  *«***♦  ►*  +  ■****•*♦* 

300  COMMON  /BODY/  R<200> ,S<200>,UB<200>,A<200) ,RI (200), NP 

320  COMMON/ 1P0T / IPOTF 

400  DIMENSION  Y  ( 1  0 ) , DY (  1  0 ) ,PRMT < 1 0) 

600  EPS  1  =  1 . OE  -  4 

700  EPS2= 1 .OE -5 

800  URITEU.I) 

900  1  FORMAT ( 1H1 ) 

1100  TEND= 100.0 

HOO  2  FORMAT ( 2 15 , F 1 0 . 5 ) 

1  720  UR  I TE ( 6 , 1 5 ) 

1725  RE AD (5,/)  NP, IPOTF .SCALE, IAN 

1727  NP 1 =NP- 1 

1730  15  FORMAT (26H  INPUT  NP , IPOTF .SCALE , IAN ) 

1740  IF( IAN.GT .0)  GO  TO  19 

1750  URITE16, 16) 

1760  16  FORMAT ( 27H  INPUT  THE  BODY  STATIONS  S) 

2000  RE  AD ( 5 , / )  (S( I ) ,  1  =  1  ,NP  ) 

2025  URI TE ( 6 , 1 7 ) 

2050  1  7  F ORMAT ( 27H  INPUT  THE  BODY  OFFSETTS  R) 

2100  R E A D ( 5 , / )  < R( I ) , 1  =  1 ,NP) 

2200  5  FORMAT ( 6F 1 0 . 5 ) 

2300  IF (SCALE. EQ. 0.0)  SCALE= 1  .0 

2400  DO  10  1=1 , NP 

2500  S( I  )  =  S( I  )/SCALE 

2600  R<  I  )  =  R(  D/SCALE 

2700  10  CONTINUE 

2950  13  CONTINUE 

2960  UR  I TE (6 , 1 8  ) 

2970  18  FORMAT  <  30H  INPUT  THE  SURFACE  VELOCITIES) 

3000  READ(5, / )  < A ( I > , I  =  1  , NP 1 > 

3100  3  FORMAT (6F13.5) 

3200  UP  < 1 >  =0 . 0 

3300  DO  9  1  =  1  ,NP 1 

3400  J«I+1 

3500  UB ( J )=A ( I) 

3600  A( I )=0.0 

3700  9  CONTINUE 

3800  CALL  ARC 

3810  GOTO  100 

3900  URI TE<  6 ,6 ) 

4000  6  FORMAT  ( 753H  R  S  LIB  A  \ 

\  R  I  ) 

4200  DO  7  1=1, NP 

4300  URITE<6,8)  R ( I) ,S ( I ) ,UB ( I ) , All) ,RI (I) 

4400  7  CONTINUE 


>5 


T 


- 


4500 

8 

FORMAT ( 5E 1 2. 4 ) 

4600 

100 

CALL  SPl INE(A,UB,S,NP) 

4700 

CALL  SPLINE ( A .RI ,R,NP) 

4750 

19 

CONTINUE 

4800 

N=2 

4850 

I  ST  EP  =  1 

4900 

PRhT ( 1 >>0.0 

49)0 

IF ( IPOTF . EQ. 0  )  GO  TO  14 

4920 

WRITE (6,20) 

4925 

20 

FORMAT (38H  PUT  IN  LOCATION  OF  FRONT  STAG.  POINT) 

4930 

READ15,/)  SSTP 

4940 

PRMT ( 1 )=SSTP 

4950 

14 

CONTINUE 

5000 

PRMT <  2)  =  TEND 

5050 

WRITE ( 6, 2 1 ) 

5100 

21 

FORMAT ( 27H  READ  IN  INITIAL  STEP  SIZE) 

5150 

READ ( 5 f / )  PRMT ( 3 ) 

5200 

PRNT ( 4  >  =  EPS  1 

5300 

PRMT ( 5 )  =EPS2 

5400 

PRMT (6)=0.97 

5500 

PRMT ( 7) =2.0 

5600 

PRMT ( 8 )  =0 . 1 

5650 

PRMT ( 1 0 )  =  1 . 0 

5700 

IF ( IPOTF .EQ . 1 )  GO  TO  11 

5800 

C=Q . 3225 

5900 

D=0 . 1 605 

6000 

GO  TO  12 

6100 

1 1 

CONTINUE 

6200 

C=0 . 445 

6300 

D=0 . 240 

6400 

12 

CONTINUE 

6500 

Y(1)=C 

6600 

Y  ( 2  )=D 

6700 

URITE16, 1 ) 

6800 

URITE(6,4) 

6900 

4 

FORMAT (43H  X, HI ,R0B,U8, C ,[if DUBDS, DU8DS2,H,CF , DELST , T / ) 

7100 

CALL  MTPC(Y,DY,PRNT,N,ISTEP, IAN) 

7200 

END 

7300 

SUBROUTINE  MTPC ( Y , DY, PRMT , N , ISTEP , I AN ) 

7400 

DIMENSION  Y< 10) , DY ( 1 0) , Z ( 1 0 ,5 ) ,PRMT( 10) 

7500 

NOUT=0 

7600 

I E  N  D  - 1 

7650 

NZ«1 

7800 

HI=PRMT (3) 

7850 

X=PRMT ( 1 ) +H I 

7900 

HO=HI 

8000 

ISTEP=1 

8500 

NOUT -NOUT ♦ 1 

8600 

18 

CONTINUE 

8710 

DO  1  1=1, N 

8720 

Z< 1 , 3  )  =  Y ( I ) 

3(> 


0735  Z  ( 1 , 1 ) =Y ( I ) 

8740  1  CONTINUE 

8750  K  =  1 

8760  14  CONTINUE 

8765  CALL  CALC ( Y ,DY ,X ,HI , N , IAN  ) 

8770  E  =  0 . 

8780  DO  2  1=1, N 

8790  Y  <  I )  =Z(  I , 3  )  +  HI  »DY  < I ) 

8800  E - A  BS 1 Z ( 1 , 1 ) -Y ( I ) ) +E 

8810  2  CONTINUE 

8812  IF(K.EQ.I)  GO  TO  17 

8815  IF < E . LT . 1 .OE-4)  GO  T03 

3820  IF(K.GT.IO)  GO  TO  4 

3825  17  CONTINUE 

8830  K=K+ 1 

8840  HO  13  1  =  1  ,N 

8850  Z (I f  I  )  =  Y  ( I ) 

8860  13  CONTINUE 

8880  GO  TO  14 

3890  4  CONTINUE 

8900  HI=HI/2 . 

8910  HQ=H I 

8912  DO  21  1=1 , N 

3914  Y( I )*Z( I , 3) 

8916  21  CONTINUE 

8920  IF ( HI . L T . 1 .OE -5 )  URITE(6,19) 

8930  19  F ORHAT ( 308  NO  CONVERGENCE  AT  FIRST  STEF) 

8’40  K= 1 

8950  GO  TO  14 

3960  3  CONTINUE 

8970  LALL  CALC( Y, DY ,X ,HI , N , IAN) 

8980  DO  20  1  =  1, N 

3990  2(1  ,1  )=Y( I ) 

7000  Z(1 ,2)=DY < I ) 

9010  20  CONTINUE 

9050  CALL  OUTPUT i X ,Z  r PRHT , I  END , HO , H I ,NOUT ,NZ ) 

11100  10  CONTINUE 

11300  C=HI7H0 

11400  CI=C*C 

11500  C2= 1 . - C 1 

11600  C3=HI*t1.+C) 

1  1700  C4=3.+2./(C*C1 ) 

11800  HI2-H1/2. 

1  1900  DO  5  1  =  1, N 

12000  Y  < I )«C1*Z(I,3)+C2*Z(1, 1 )*C3*Z(1 ,2) 

12100  5  CONTINUE 

12300  X=X+Hl 

12400  CALL  CALC(Y,DY,X,HI,N,1AN) 

12500  E*0.0 

12600  DO  6  1=1, N 

12700  Z< 1 ,5  )=Y ( I ) 


12800 

Y 1 1  )  =  Z  1 1 , 1  )+HI2*(DY(I)+Z(It2>) 

1 2900 

AE«A8S1Z<I,5>-Y<I)>/C4 

13000 

IF(AE.GT.E)  E-AE 

13100 

6 

CONTINUE 

1  3300 

IF ( E . GT . PRHT ( 4 ) >  GO  TO  7 

13400 

CALL  CALC(Y,DY,X,HI,N,IAN) 

1  3500 

DO  9  I  =  1 ,  N 

13600 

2 ( 1 , 3 )=Z ( 1 , 1  ) 

13700 

*— » 

M 

6— « 

fSJ 

13800 

Z ( 1 , 1 )  =Y( I ) 

13900 

Z <  1 ,2) =DY ( I ) 

1  4000 

9 

CONTINUE 

14100 

I  STEP=2 

1  4200 

N0UT=N0UT+1 

1  4300 

CALL  OUTPUT ( X,Z,PRNT, IEND  ,HO,H I , NOUT ,NZ > 

14400 

IF ( IEND.EQ.-1 )  GO  TO  15 

1  4500 

IF1E.LT. PRNT<5>)  GO  TO  8 

14600 

H0  =  H  I 

1  4700 

GO  TO  10 

1  4800 

T 

/ 

CONTINUE 

15000 

IF ( ISTEP.EQ. 1 )  GO  TO  1 1 

15100 

A= 1 PRNT (6) *PRHT  <  4 )/E )**( 1 ./3. ) 

15200 

X  =  X  -  H I 

15300 

H1=A*HI 

15400 

IF1HI.LT. 1.0E-7)  GO  TO  12 

15500 

GO  TO  10 

15600 

11 

CONTINUE 

15800 

X=X-2.*HI 

15900 

HI  =HI/2 . 

16000 

IF1HI.LT. 1 .0E-10)  GO  TO  12 

16100 

1F1E.LT. I.0E-10)  E=PRHT 1 5 ) 

1  6200 

8 

CONTINUE 

1  6400 

A  = 1 PRHT ( 5 ) +PRHT (7)/E)**(1./3.) 

16500 

HO  =  HI 

16600 

HI=A*HI 

16700 

IF1HI.GT.PRNT18))  HI=PRHT ( 8 ) 

16800 

GO  TO  10 

16900 

12 

CONTINUE 

16920 

PRNT (10)=-) .0 

1  6940 

PRNT1 1 >=0.0 

1  6960 

CALL  OUTPUT(X,Z, PRNT, IEND, HO, HI, NOUT ,NZ) 

17000 

UR I TE 1 6 , 1 6 )  X .HI 

17100  16 

FORMAT 1/27H  STEP  SIZE  TOO  SMALL  X=,E12.4,4H  HI = ,E 

17200 

RETURN 

17300 

15 

CONTINUE 

1  7400 

RETURN 

17500 

END 

17600 

SUBROUTINE  ARC 

17700 

COMMON  /BODY/  R (200) ,S< 200 ) ,UB( 200) , A< 200 > ,RIl 200) ,NP 

17800 

RI11)*0.0 

17900 

A 1 1  ) =0 . 0 
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18000 
18100 
1  8200 
18300 
18400 
18500 
18600 
1  8700 
18800 
18900 
19000 
19100 
19200 
19300 
19400 
19500 
19600 
19700 
1  9800 
19900 
20000 
20100 
20200 
20300 
20400 
20500 
20600 
20700 
20800 
20900 
21000 
21100 
21200 
21300 
21  400 
21500 
21600 
21700 
21800 
21900 
22000 
22100 
22200 
22300 
22400 
22500 
22600 
22700 
22800 
22900 
23000 


IFLAT=2 
IFLAT 1 =  1 
DO  1  1=2 ,NP 

IF (S ( I ) -S ( I- 1  I.EQ.O.O)  GO  TO  1 
I F  L  A  T  =  I 
GO  TO  2 

1  CONTINUE 

2  CONTINUE 

IFUFLAT.LT. 3)  GO  TO  3 
IFLAT 1 =1FLAT-1 
DO  4  I =2, IFLAT 1 
A ( I )=0 .5* (R( I ) +R( 1-1 ) ) 

RI ( I )=A( I ) 

4  CONTINUE 

3  CONTINUE 

X 1 1 =R ( IFLA  T 1  ) 

DO  5  I = IFLAT , NP 

F  =  (SU)-SU-1)>*SQRT<1  .  +  ( (R(I)-R(I-I ) >/ <S< I )-S< 1-1 > > >**2> 
RI ( I ) =0. 5* ( R ( I ) +R ( I  - 1 ) ) 

X I =X  1 1 *F 

A ( I )=0 .5* (XI+XI 1  ) 

XI 1 =XI 

5  CONTINUE 
RETURN 
END 

SUBROUTINE  SPL I NE ( X f Y ,C , N ) 

DIMENSION  X(200),Y<200), Ai200) ,B(200) , D < 200 > , C < 200 ) 

N 1 =N- 1 
N2=N-2 
DO  1  1=1 ,N2 
J=IM 

HI  =  X( J)-X( J-l  ) 

H1 1 =X( J  +  l )  -XiJ) 

D(I)  =  3.*(((Y(J+1 ) -Y<  J ) )/H1 1 ) - <  Y  L  J )  -Y(J-1 ) » / HI )  /(HI+HI1  ) 
B<I)=HI1/(HI+HI1 ) 

A(I)  =  t  .-BUI 
A( I)*0.5»A(I) 

B(I)=0.5»BU> 

1  CONTINUE 
Ci  1  )  =  1 .0 

B( 1 )=B( 1 )/C(  1  ) 

DO  2  I =2, N2 
C ( I )  =  1 .0-A ( I ) *B ( I- 1 ) 

B ( I )  =B (I )/C( I ) 

2  CONTINUE 

A  < 1 ) =D  < 1  )/C«1  ) 

DO  3  1=2, N2 

A  ( I  )  =  (  D  ( I )  A(  I )  ♦  A  ( I  -1  >)/CU> 

3  CONTINUE 
N3=N2  -1 

DO  4  1  =  1  ,N3 


23100 
21200 
2  3300 
23*00 
23500 
23600 
23'00 
23800 
23900 
2*000 
2*100 
2*200 
2*300 
2**00 
2*500 
2*550 
2*600 
2*200 
2*800 
2*900 
25000 
25100 
25200 
25300 
25*00 
25500 
25&00 
25200 
25800 
25900 
26000 
26100 

26200 

26300 

26*00 

26500 
2651  0 
26520 
26522 
2652* 
26526 
26528 

26530 

26531 

26532 
2653* 
26536 
26538 


J=N2-I 

C(J>»A(j)-Bu>«cu+n 

*  CONTINUE 
CiNI )*0.0 
C<N)*0.0 
DO  5  1=1, N3 
J=N2- 1 
K=JU 
C  < K  >  =C( J > 

5  CONTINUE 
YU)  =0.0 
RETURN 
END 

SUBROUTINE  POTF l X ,U ,BU  ,DU2 , RO ,PRO, 1  AN ) 

COHNON  /BODY/  R(200),S(200),UB(200),A(200>,RI(200),NP 
IRIAN. GT.O)  GO  TO  3 
N  =  NF 

DO  1  1=1, NP 
IF  vX.GE.A1 1) )  GO  TO  1 
N*  I 

GO  TO  2 

1  CON  T I  HUE 

2  CONTINUE 

H J  =  A i N 1  - Ai N- 1 1 
X3  =  <  HJ»  *2  )/6 . 

X4=(AiN)-X)/HJ 

X5MX-ATN-1 ) )/HJ 

Y1=( <AiN)-X)**2)/(2.*HJ) 

Y2=((X-A(N-1))**2)/(2.*HJ) 

XI  * ( A ( N ) -X )*Y 1 /3. 

X2=(X-A\N-I ) )*Y2/3. 

U*S»N-I  )*X1*S(N)*X2*(UB<N-1 )-S<N-1 )*X3)*X4MUB(N)-S(N)*X3)\ 

\  •  X  5 

R  0  =  R  <  N -  1 )«X1*R(N)*X2+(RI(N  1  )-RlNU )*X3)*X*  +  (RI(N)-R(N)*X3)\ 
\  *X5 

D0»-S(N-1 )*Y1 ♦  S(N>»Y2+(UBiN>  -UB<N-1 ))/HJ-(S(N)  -Si N  -1 ))*HJ/6. 
DRO= -R(N-1 >*YRR<N)*Y2MR1(N)-RI<N- 1 ) )/HJ- (R(N)-R(N-I ) ) *HJ\ 
\/6 . 

DU2  =  S i N  - 1 ) »X  *  *S i N ) *X5 
3  CONTINUE 

IFUAN.GT.  1 )  GO  TO  * 

RO=SIN(X ) 

DRO=COS(X) 

U=t .5*R0 
DU=1 .5«DR0 
DU2=-U 
RETURN 
4  CONTINUE 
R0=1  . 

DR0=0 . 

U=2.«SINiX) 


26540  DU=2 . *C0S 1 X ) 

26542  DU2=-U 

26600  RETURN 

26200  END 

26800  SUBROUTINE  OUTPUT ( X , Z ,PRMT , IEND ,H0 , HI , NOUT ,NZ ) 

26900  DIMENSION  Z( 10,5)  ,PRNT ( 10) 

27000  COMMON  /BL1/  J,SL ,SM ,H , DST AR , HE ,HT ,UB ,DUBDS  ,R0B , DROBDSY 

\ , DUBDS2 

27)00  COMMON  /BL 2/  Y < 4 > .RATIO, G, DTPS 

27200  REAL  J 

27225  1F1N0UT.EQ.1  )  Y ( 1 )  =  1 . 

27250  NZ 1 =500 

27300  Y( 4 )*DTDS 

27400  Y  <  2  >  —  Y  < 1 ) 

27500  Y ( 3  >  =HI 

27600  YU)  =  -RATIO 

27820  XRAT 10= -RAT  10 

27830  IFURATIO.LT. 0)GOTO  10 

27840  T=SQRT (XRAT 10) 

28000  CF=2. *U6*SL/T 

28100  DELST=H*T 

28150  IF ( PRMT 1 101.LT.0.0)  GO  TO  12 

28200  IF ( CF . GT . 0 . 1  )  GO  TO  12 

28275  PRMT ( 3 ) =PRMT 131/10. 

28280  PRMT (10)=-'  .0 

28290  12  CONTINUE 

28300  A=X-PRNT < I ) 

28400  IF1A.LT. PRNT131 )  GO  TO  1 1 

28500  PRMT ( 1  )  *X 

28900  2  CONTINUE 
28950  DUB2-DUBDS2 

29000  URITE<6,1>  X ,HI ,ROB ,UB ,Z ( 1 , 1 ) , Z (2, 1 ) .DUBDS.DUB2  ,H,CF\ 

N.DELST.T 

29100  1  FORMAT ( 1 2E 1 0 . 3 ) 

29150  11  CONTINUE 

29200  8  IF1CF.LT. 0.)  IEND=-1 

29300  IF ( X . GT . PRM T 1 2  ) )  IENP=-1 

29400  RETURN 

29450  10  URITE16.201X 

29470  20  F ORMAT ( 1 X ,//, 3X , 1 OHCF  NEGATIV  ,3X,3HX=  .F10.4) 

29480  STOP 

29500  END 

29600  SUBROUTINE  CALC  1 Y , DY . X ,BX ,N , IAN ) 

29650  DIMENSION  YINI.DYiN' 

29700  COMMON  /BL 1 /  J. SL , Srt.H, DS TAR ,HE , H T ,UB, DUBDS ,ROB, PROBDSX 

\ ,  D  U  B  D  S  2 

29750  COMMON  / BL 2 /  Z 1 4 ) .RAT  10, G ,DTDS 

29900  REAL  J.L.M.MPD.MPC 

29950  C  =  Y  ( 1  ) 

29960  D  =  Y ( 2  > 

30000  CALL  P0TF1X.UB, DUBDS, DUBDS2, ROB, DROBDS, IAN) 
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30100 

30200 

30300 

30400 

30500 

30600 
30700 
30800 
30900 
31000 
31 100 

31200 

31300 

31400 

31500 

31600 

31700 

31800 

31900 

32000 

32100 

32200 

32300 

32400 

32500 

32600 

32700 

32800 

32900 

33000 

33100 

33500  10 

33600 

33700 

33800 

33900 

34000 

34100 

34200 

34300 

34400 

34500 

34550 

34600 

34660 

34700 

34800 

34900 


J*.  32832+(  -.  31 897+.  1 4931  *0*C  +  P*  (.  18760-.  14947*0 
1*1. 747+3. 5*C 

ESC*.  20069* (-. 054971  +(-.18242+1. 17580+1. 04J227+(-. 055451 +\ 
0095503* 

»0*0*C)*C)*0*C 

E5P= .060072+ 1 .26478+( - . 29465+ 1 -. 1 1 6 1 1 ♦(. 1 5200- .029892*0 *C)  \ 

\*0*0* 

$C 

ESP2*-.  1 1365+1  .089165+1  .1  221  5  +  (  -.13673+. 031  21  9  +0  *0  * C  )  * C 
ESD3*. 020669+ (-.049679+ (.0401 41 -.010879*0  +C1+C 
E=ESC+(ESD+(ESD2+ESP3*D)+B) *  D 
EPD=ESD+ 1 2 . +ESD2+3 . +ESD3+P)  *P 

EPCC  = -.054971+ (-.36484+1 .52740+1. 17291 +1 -.27726+. 057302*0 \ 
\*C)*C)*C 
».'*C 

EPCD*.  26478 +  (-.  58930+ (-.  34833+ 1.60800-.1 4946*0  *0*0*C 
EPC D2=.089 165+1 .24 43+1 -.41019+. 12488*0 *C)*C 
EPCD3*-. 049679 +(.080282*. 032637+0 *C 
EPC  =  EPCC+ 1 EPCD  +  l E  PC  D2+EPC  D3*  D ) *  D ) *0 

GSC=. 12733+1 -.050050+1*. 10089+1  .  1391 5-.  035073  :*C)*C>*0  +  C 

GSP=.  04905  4  +  1 . 13403  +  1  -.  23383+ .072643  *0*0  +  C 

GSD2*-. 059455+ 1.09391 6-. 037702*C)*C 

G=GSC+(GSD+GSP2*D)+D 

GPD=GSD+2.*D+GSB2 

G P C C  =  -.050050+ (-.201 78  +  1 .41 745-. 14029*0  *C) *C 

GPCD =.13403+1*. 46766+. 21793*0  *C 

GPCD2= .09391 6*. 075404  *C 

GPC  =  GPCC  +  (GPCP+GF'CD2*D)*D 

8= -30 . *0 

flPC  =  0. 

«PP=-3o. 

HE=E/G 

G2=G*G 

Srt=G2*M 

SL=G*L 

H=J/G 

Al  =  1 .3713+1 1 .2124+1 1 .8754+1-1 .8268 + .60524*0 *C ) *C)*C 
A I=AI+D*< -.83784+1*1 .9728+13.3493-1  .2919*0*0*0 
A1=AI+B*D*< .94797+1-1  .5608  +.69570*0*0 
DSTAR=G*AI 
GH=G+N 

HEPD=(EPD*G-GPD *E)/G2 
HEPC  =  (EPC  +G-GPC  *E  > /G2 
SNPP  =  2.  *G«  *GPD  +  G2*f1PD 
SNPC=2.+GM*GPC 
A02D=0. 

IFIABSIDUBDS) .LE.1 .0E-2)  GO  105 
IF(IPOTF.EO.O)An2Ii=SH*UB*PROBDS/(DUBDS*ROB) 

DTDS=2.*( SL+SM* 1 2.+H) +A02D1/UB 

RATIO=SH/DUBBS 

GO  TO  6 


35000  5 

35050  C 

35060 

35070 

35080 

35100 

35200 

35300 

35350 

35400 

35500 

35600 

35650 

35700 

35895  210 

36000 

36050 

36100 

36200 

36300 

36400 

H 


CONTINUE 

PRIVIOUS  VALUE  IN  Z(4) 

TSIH  =  Z< 1  ) 

TSI=Z<2) 

DT  DS I =Z ( 4 ) 

TS=TSIN+(2.*<TSI-TSIN)/Z(3)-DTDSI)*<Z<3)+DX) 
TS=TS+<BTDSI/Z(3)-<TSI-TSIN)/(Z<3)*Z<3)))*<DX+Z<3) >**2 
RAT I0=-TS 

IF(IP0TF.EQ.0)A02D=-TS*UB+DR0BDS/R0B 

DTDS=2.+(Sl+SM*(2.+H)+A02D)/UB 

CONTINUE 

R1  =-2.:*BUBBS+  ( SL+SN+( 2.  +H) ) /UB 
IF(IF'OTF.EQ.O)AQ2D=SN*IiROB[iS/ROB 
R1=R1-2.*A02D+RATI0*BUBBS2 
R2=HE  »( SL  +  SN*  < H - 1  > ) 

DS=SMPC*HEPB-SMPB*HEPC 

R2=- <2 .+DST AR  R2 )/( RATIO -»UB  ) 

DY(  1  )  =  (R1  *HEPD-R2;*SMPD)/DS 
HY ( 2 )  =  ( R2*S«PC  R 1  +HEPO/DS 
RETURN 
END 
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AITOL/FYS,  J.  Olsen 

Nat  tonal  Hut  o.i u  ol  Standards 
I  1’ .  S .  K 1 e Hanoi  I 

MARA!' 

1  lUv  ot  Ship  R&D 
l  Library 

NASA,  IIQS/Llb 

NSF/F.ng.  Lib 

LC/SCl  &  TECH 

nor/i.fb.  tad-A‘M.1 

MMA 

1  Nat’ l  Marti  Into 
Res.  Con, 

1  Library 

La  I  .  Inst.  Techno  1  ogv 
l  A.  Costa 

('•it  I  .  lust.  Techno  logy 
•  lot  Propti  1 a  ton  Irfih.  , 
Pasadena,  Citl. 

1  L.M.  Mack 

Case  Western  Hn  I  vet  s  I  t  v 
1  K,  Reshotko 

Uttlv.  ot  Cit  l  l  t  Of  U  lit  ,  Ik'S 
Auge  l  es 
I  S.J.  Marker 
I  A.R.  Wait /.an 


l  Illinois  Inst,  of  Tech. 

I  M.V.  Morkov Irt 

l  tin  tv.  ot  Iowa 

l  L.  tandweber 

l  Venn  State  University 

I  J..1.  Elsenhuth 

.'  Penn  State  Unlv./AKL 

1  B.  Parkin 
l  C.  Lauchlc 

I  l’r  l nc ot  on  Un  1  vers  1 1  v 

l  S . I .  Cheng 

Massachusetts  Inst,  ot  fecit 
l  S.A.  Or  sr.ag 
1  P,  Loehoy 

l  I'nlv.  ot  Minn.,  SAKAlll. 

I  R.  Arndt 

I  Virginia  Polytechnic  Inst, 

and  Stale  Unlv. 

1  W.S.  Sarto. 

1  H.P.  Tel  touts 

l  tlnlv.  ot  Rhode  l  si. -nut 

I  K.M.  White 

l  SIT  Uav  idson  lab., 

1  .1 .  liras  I  1  n 

l  Ailv.  Tech.  Center,  Inc., 

Texas 

I  C . S .  We l 1 s 

I  MOW  1  ndiist  i  I  i'm  ,  I  nc  . 

I  O.C.  Wilcox 

1  Lockheed  Alici.it  t,  Sunnvvab 

I  R .  Wit  i  d 

.’  Mclkinnelt  llenglas  Corp. 

1  Lib 
l  T.  Cebecl 


Cop  i  os 

Cop las 

Code 

Name 

l 

Physical  Dynamics,  Inc,, 

1 

1902 

C.  Maidanlk 

Torrance,  Cal. 
t  W.  Kalgh 

1 

1903 

U.  Cher lock 

1 

1905 

K.K.  Blancardl 

•) 

Rand  Corp.,  Santa  Monica, 

1 

192 

Cal. 

1  C.  Caz ley 

1 

191 

R.C.  Olson 

l  J.  Aroesty 

1 

1933 

J.P.  Lee 

3 

Rockwell  Intern. 

1 

1942 

.1 .  T.  Shen 

Au t one t lcs  Croup 

1 

1942 

R.E.  Bowers 

1  D.B.  Moody 

1 

1942 

M.J.  Casarrlla 

l  C.  Jennings 
l  B.  Ca  rnt  Ichael 

l 

1946 

.l.K.  Spina 

1 

West  lug house  Klee  trie  Corp. 

10 

5211.1 

Report  s  lllst  r  Unit  Ion 

Annapolis,  Mil. 

1 

522.1 

line  lass  1 1  led  Lib  (Cl 

1  l.lbrary 

1 

522.2 

line  lass  ll  led  Lib  (A) 

1  West Inghouse  Research  Co r p . 
Pittsburgh,  Pa. 


l  F, 

.  R.  Co  hi  so  hm  led 

CK.NTKK 

DISTRIBUTION 

Cop l es 

Code 

Name 

l 

15 

t 

1504 

V.J.  Monaco l la 

1 

1522 

J.H.  Robinson 

l 

1 524 

l 

1532 

C .  Do  ha  y 

1 

l  54 

W.B.  Morgan 

l 

1  552 

J.H.  McCarthy 

l 

1552 

N.  droves 

1 

1552 

T.  T.  Ilnang 

30 

1  552 

C.  von  kei o.rok 

l 

184 

J.W.  Sc hot 

l 

184  3 

C.  Dawson 

l 

19 

M.M.  Sevtk 

1 

1901 

M.  Strasberg 

DTNSRDC  ISSUES  THREE  TYPES  OF  REPORTS 

1.  DTNSRDC  REPORTS,  A  FORMAL  SERIES,  CONTAIN  INFORMATION  OF  PERMANENT  TECH¬ 
NICAL  VALUE.  THEY  CARRY  A  CONSECUTIVE  NUMERICAL  IDENTIFICATION  REGARDLESS  OF 
THEIR  CLASSIFICATION  OR  THE  ORIGINATING  DEPARTMENT. 

2.  DEPARTMENTAL  REPORTS.  A  SEMIFORMAL  SERIES.  CONTAIN  INFORMATION  OF  A  PRELIM 
INARY,  TEMPORARY.  OR  PROPRIETARY  NATURE  OR  OF  LIMITED  INTEREST  OR  SIGNIFICANCE. 
THEY  CARRY  A  DEPARTMENTAL  ALPHANUMERICAL  IDENTIFICATION. 

3.  TECHNICAL  MEMORANDA,  AN  INFORMAL  SERIES.  CONTAIN  TECHNICAL  DOCUMENTATION 
OF  LIMITED  USE  AND  INTEREST.  THEY  ARE  PRIMARILY  WORKING  PAPERS  INTENDED  FOR  IN 
TERNAL  USE.  THEY  CARRY  AN  IDENTIFYING  NUMBER  WHICH  INDICATES  THEIR  TYPE  AND  THE 
NUMERICAL  CODE  OF  THE  ORIGINATING  DEPARTMENT.  ANY  DISTRIBUTION  OUTSIDE  DTNSRDC 
MUST  BE  APPROVED  BY  THE  HEAD  OF  THE  ORIGINATING  DEPARTMENT  ON  A  CASE-BY-CASE 
BASIS 


1  0,M.  Griffin 

NUSC,  N'LONLAB  j 

1  R.A.  Skop 

1  rl.  ScMoemer  j 

RDA 

1 

NSWC,  Dahlgren/ Lib 

NA 

2 

NSWC,  White  Oak 

1  Dr.  B.  Johnson 

1  Dawson 

1  Nav.  Sys.  Engr.  Dept 

1  Tech  Library 

1  Library 

VPGSCOL 

1 

NAVSEC,  N0RVA/660.03  Blount 

1  Library 

X  T.  Sarpkaya 

1 

NCEL/Code  131 

1  J.  Miller 

1 

NAVFAC/ Code  032C 

DC 

1 

NAVSHIPYD  PTSMH/Lib 

SC/712,  D.  Humphreys 

1 

NAVSHIPYD  PHILA/Lib 

VSEA 

1  SEA  03 

1 

NAVSHIPYD  NORVA/Lib 

2  SEA  03D 

1 

NAVSHIPYD  CHASN/Lib 

2  SEA  05H,  A.R.  Paladino 

1  SEA  05R 

2  SEA  312 

1 

NAVSHIPYD  LBEACH/Lib 

47 
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